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We use a partially Gutzwiller projected BCS d-wave wavefunction with an antiferromagentic 
weighting factor to study the ground state phase diagram of a half filled Hubbard-Heisenberg model 
in a square lattice with nearest neighbor hopping t and a diagonal hopping t' . The calculations 
are carried out by using variational Monte Carlo method which treats the Gutzwiller projection 
explicitly. At large on-site Coulomb interaction U , the ground state is antiferromagnetic. As U 
decreases, the ground state becomes superconducting and eventually metallic. The phase diagram 
is obtained by extensive calculations. As compared to the strong effect oiU/t, the phase boundaries 
turn out to be less sensitive to t' /t. The result is consistent with the phase diagram in layered organic 
conductors, and is compared to the earlier mean field result based on the Gutzwiller approximation. 
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I. INTRODUCTION 

\ High-temperature superconductivity remains to be an 
exciting and rich field. One of the interesting proposals 
is Anderson's resonating valence bond (RVB) stateJ^i^ 
In the RVB theory, the parent compound is an insulator 
at half electron filling, or one electron per Cu-site, and 
chemical doping is essential to introduce charge carriers 
to lead to superconductivity. The mathematics of the 

\ RVB theories therefore is in a Hilbert space which com- 
pletely projects out the on-site double-occupied electron 
states. At the half filled, there is exactly one electron per 
lattice site, and the charge degree of freedom is totally 
frozen, resulting in a Mott-insulator. 

Another interesting class of materials in the context 
of strongly correlated systems is the layered organic 
conductors ^i^i^ii which may undergo a phase transition 
from an insulator to a superconducting (SC) state by 

\ applying pressure Since these materials are effectively 
at half fillingji^, the phase transition is due to the com- 
petition between the Coulomb interaction and kinetic 
band width, the latter of which is tuned by pressure in- 
stead of chemical doping. There have been several re- 
lated theoretical works on layered organic superconduc- 

■ tors in recent years ^^^^^i^^ii^^ii^ia The math- 

\ ematics of the SC state may be described by a partially 
Gutzwiller projected BCS state,— instead of the com- 
plete projection as in the RVB theory. We shall re- 
fer to this partially Gutzwiller projected BCS state as 
a Gossamer superconductor, a phrase first introduced 
by Laughlin^i^ originally in the context of high tem- 
perature superconductors. Gossamer superconductivity 
refers to those SC states with a dilute superfluid-density. 
The partial Gutzwiller projection allows charge fluctua- 
tions even at half filling. One of us^" proposed that in 
this case an effective model is the Hubbard-Heisenberg 
model which includes the standard kinetic energy, the on- 
site Coulomb repulsion, as well as the anti-ferromagnetic 
spin-exchange. The idea was applied to the study of k- 



(BEDT-TTF)2X by Can et al^, where the Gutzwiller 
approximation was used to replace the partial Gutzwiller 
projection by a set of renormalized factors and the re- 
sulted renormalized Hamiltonian was then studied by a 
mean field theory. The finding is a phase-diagram dis- 
tinguishing three phases: normal metal, superconductor 
and anti-ferromagnet. 

In this paper, we shall study the phase diagram of an 
effective Hubbard-Heisenberg model by using variational 
Monte Carlo (VMC) method. The order-parameters 
for the d-wave superconductivity and for the anti- 
ferromagnetism are calculated directly. We obtain a 
phase diagram consistent with the experiments, provid- 
ing further support to the scenario of the Gossamer su- 
perconductivity to describe the layered organic conduc- 
tors. Interestingly, our numerical calculation of the SC 
order parameter suggests a relatively high superfluid den- 
sity near the phase boundary to the AFM insulator. The 
results from our VMC calculations also provide support 
to the earlier mean field results based on the Gutzwiller 
approximation,"'^^ although we find a less sensitive role of 
t'. 



II. MODEL, TRIAL- WAVE-FUNCTION AND 
METHOD 

We study a Hubbard-Heisenberg model in a 2- 
dimensional lattice illustrated in Fig. [TJ The Hamil- 
tonian is given by 

= ^ Un^^Uii ~ ^ tijcl^Cja + h.c. + '^JSi- Sj{l) 
* (hi) 

Here Cj^ is the electron annihilation operator of an elec- 
tron with spin a on site i, Si is the spin- 1/2 operator 
at site i, and rii^ — cl^Cia- The sum is over the 

nearest neighbors (n.n.) pairs on the square lattice, and 
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FIG. 1: (Color online) Illustration of the lattice for Hamil- 
tonian (1) studied in this paper, t and J are hopping and 
spin exchange coupling between nearest neighbor pairs (solid 
lines), and t' is the hopping integral along a diagonal direction 
(dashed lines). 



the sum over is over both the n.n. pairs and the 

diagonal bonds (dashed hnes in Fig. 1). We set the n. n. 
hopping tij=t = l as the energy unit, and fix the spin 
exchange J as J/t = 0.5, and treat the diagonal hopping 
integral tij = t' and the on-site repulsion U as tuning 
parameters. In our numerical calculations, we consider 
a L by L lattice and use a periodic boundary condition 
along the x-direction and antiperiodic boundary condi- 
tion along the y-direction. 
Our trial wave-function reads 



I*' 



SfS' 



(2) 



where \'^n) is the BCS-wave function projected to the 
subspace with the fixed number of particles as defined 
in Eq. (3) below. In Eq. (2), we introduce two variational 
parameters a and /?, to control the partial Gutzwiller 
projection and the AFM correlation, respectively. The 
BCS-wave function in the fixed particle formalism has 
the following form in real space: 



where Rja- is the spatial position of an electron with spin 
a at the lattice site j, and the sum is over all the pairs of 
a spin-up electron at the site j and a spin-down electron 
at site /. Here a(r) is the amplitude of the wave function, 
which is the Fourier transform of a{k) = Vk/uk, with 
and Vk given in the usual BCS wavefunction: 



I* 



BCS) 



n( 



)|0) 



(4) 



We haveii^S; 



a{r) — flfc cos(fcr) 
fc 

A(fc) 



(5) 



Following the previous literature on the pairing sym- 
metry for the modeliSiiii, we focus here on the d^^-y^- 
wave pairing state, where A(fc) and ^fc have the following 
forms, 



A(fe) = A(cos(/cj;) - cos{ky)) (6) 
-2t(cos(fc^) -I- cos(fcj^)) — 2t'^ cos{kx + ky) — /i (7) 



where A, jj, and are variational parameters in the the- 
ory. Note that t'^ ^ t' in general due to the spin coupling 
term in the Hamiltonian. The advantage of the above 
trial wave function is that the SC and AFM order can 
be treated on an equal footing. It turns out that a small 
value of P improves the energy of the SC state, while 
a sufficiently large value of (3 leads to AFM long range 
ordering. We measure the staggered magnetization to 
quantitatively study the the AFM phase. 



m 



\ 



N{N - 



(8) 



where TV is the number of the lattice sites, and the sum 
is over all the N{N — 1) pairs between sites i and i -f r 
on the lattice. To measure the SC long range order, we 
introduce a pair correlation function: 



(9) 



where hi^i+r is a spin singlet bond between the two sites 



|0) (3) ^ 



and 



T, and r — ±x, ±y, Ty 



for 



zbi, and 



Ty = ±1 for T = ±y. Fi describes a d-wave singlet bond 
around the site i. The off-diagonal long range order pa- 
rameter for the d-wave pairing can be measured by the 
quantity at i? — > oo, 



(10) 



where the sum is over all the lattice sites. In our cal- 
culations on the finite size systems, we choose K ~ 
{L/2,L/2), the largest displacement on the lattice of L 
by L with L upto 10. 
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FIG. 2: (Color online) Variational ground state energy Eg of 
Hamiltonian (1) as a function of U jt for J/t = 0.5 and various 
values of t' jt. The inset is an enlarged figure for the energy. 



To simplify the variational procedure, we will fix /i 
in the calculations with the reasons given below. It 
has been argued^'^ that A and /i are not independent 
in the variational calculations for the t — J model. We 
found that the results are essentially insensitive to /i 
for the present model. In contrast, t'^ is an important 
variational parameter here. The ground state energies 
and the ground state phase is sensitive to t'^ over a 
wide parameter-range. By fixing /it, we have then four 
variational parameters (A,t^,a,/3) in our calculations 
to determine the phase diagram in the parameter space 
of U and t' . 



There are two sources of error bars in our numerical 
calculations within the variational approach. One is from 
the statistical errors, and the other is due to the discrete- 
ness of the variational parameters in our calculations. 
In our simulation, we start with several different initial 
configurations and then average our numerical measure- 
ments over those simulations. The error-bars obtained 
in these averages are found to be one order of magnitude 
smaller than the error-bars described below. We consider 
the possible values for the variational parameters and di- 
vide them into small slices. Then we perform VMC for all 
combinations of this " mesh" . After obtaining an optimal 
set of variational parameters in this mesh for a particular 
set oiU/t and t' /t, we develop a local mesh for nearby 
values of the tuning parameters. From the spacing of our 
mesh, we obtain the error-bars for the variational param- 
eters. The results and the error-bars we present in this 
paper are essentially due to the finite elements we choose 
in the variational parameters. 



III. RESULTS AND DISCUSSION 

In this section, we present our results on the varia- 
tional ground state, the corresponding variational pa- 
rameters, the SC and AFM long range orders of the 
Hamiltonian (1) in the parameter space of U/t and t' /t. 
Since the phase of the ground state is much more sen- 
sitive to the on-site repulsion U than to the diagonal 
hopping i', we focus our study on three values of t' , with 
t' /t = 0.1, 0.5, 0.9. In Fig. Hwe plot the obtained ground 
state energies as functions of U for different values of t' . 
The corresponding optimized variational parameters and 
long range SC and AFM order parameters as functions 
of U for three sets of values of t' are plotted in Fig. [31 
The simulations are carried out on lattice sizes of L=6, 8, 
and 10, as indicated in the figure. Before we discuss the 
results, we note that the spin coupling term in Eq. (1) 
is to account the virtual hopping process in the Hubbard 
model, which is derived at the large U limit. The present 
study may be of relevance to the Hubbard model only at 
large C/, but not at small U . Our main interest will be at 
large or intermediate values of [/, and the interpretation 
of our results at small U to the Hubbard model should 
be cautious. 

Before we discuss general features, we briefly discuss 
the obtained variational parameter t'^, which is to opti- 
mize the kinetic energy due to the presence of the diago- 
nal hopping integral t' . t'^ increases as t' increases, but t'^ 
is significantly smaller than t' as we can see from the first 
row in Fig. [31 At large U, t'^ becomes zero or very tiny. 
This may be understood as a result of the AFM ground 
state with commensurate wave vector (tt, tt), since a finite 
t'^ does not match the AFM state and is not preferred. 

As U increases from zero, the projection parameter 
a increases from around 0.05, indicating a graduate in- 
crease in Gutzwiller projection, while the weighting fac- 
tor parameter (3 changes little at small U, but changes 
rapidly around C/ ~ 5. The mean field pairing amplitude 
parameter A changes slowly at small [/, but increases 
rapidly starting from around U — 2, then reaches a max- 
imum at around U = 4.5 and drops at larger U. The 
ground state properties are best seen in the measurement 
of the SC order parameter cj) and AFM order parameter 
m. Qualitatively there are three regions as U increases. 
At small U, both </) and m are very tiny or essentially 
zero, indicating a metallic state. At intermediate U, (f> 
increases monotonically with U, while m remains tiny, 
indicating a SC state. As U further increases, drops 
sharply, while m increases rapidly. We may identify this 
phase as the AFM phase without the SC order. The tiny 
but non-zero values of (j) and m in the non-ordered states 
may be explained as finite size effect, although a system- 
atic scaling analyses is difficult due to the small sizes we 
have studied. The above features are qualitatively sim- 
ilar for t' /t — 0.1, 0.5, 0.9. This is somewhat different 
from the early analytic calculations by using Gutzwiller 
approximations on the projected wavefunctions, where t' 
is found to suppress the AFM phase. We note that while 
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FIG. 3; (Color online) The ground state variational parameters t'^ (Eq. [7]), a (partial projection, Eq. [2]), /? (AFM weighting, 
Eq. [H) and A (pairing amplitude, Eq. [ejl as functions of U/t for t' /t = 0.1 (left), t' /t = 0.5 (mid) and t' /t = 0.9 (right). 
J/f = 0.5 is fixed. Also plotted are the measured d-wave SC order parameter (Eq. |9} and the staggered magnetization m 
(Eq. [S]). The lattice size is L x L, with L = 6, 8, 10. The selected error bars shown are typical, due to the finite parameter 
spacing in our calculations. 



the onset for the SC-phase is similar for different t' , the 
magnitude of the SC-order parameter is much bigger for 
the t' /t = 0.5 and t' /t = 0.9. As we can see from the fig- 
ure, the largest SC order parameter (j) is found near the 
boundary to the AFM-phase. At t' ~ 0, we expect the 
model (1) to have instability towards a commensurate 
AFM state for any finite U. 

In Fig. m we plot the phase diagram of model (1) 
obtained within our variational wavefunctions. While 
the phase-boundary between SC and AFM can be found 
easily by considering one point clearly belonging to the 
AFM and one point clearly belonging to the SC phase, 
between SC and metallic phase we have to use an arbi- 



trary value to define the phase-boundary, as the onset 
of the SC-order parameter (j) is not so sharp. We choose 
4> < 0.004 as our criteria classifying the phase to be SC. 
The error-bars in this diagram reflect within which area 
we have uncertainty that a point would be in either of 
the two phases considered. Comparing the phase dia- 
gram obtained in the VMC method with the previous 
result by using renormalized mean field theoryjii they 
qualitatively agree with each other in the sense both 
give the three phases, and overall features are similar. 
However, there are two differences. First, while both of 
the methods give the transition point between the SC- 
and AFM-phases to increase when t' /t increases, in the 
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FIG. 4; (Color online) (a) Phase diagram of the ground state 
of the Hamiltonian (1) obtained from our VMC calculations: 
Metal (metallic phase), d-wave SC (superconducting phase), 
and AFM (anti- ferromagnetic phase). The region with large 
SC order parameter is indicated by a thick line and marked 
with (fimax ■ The values of the order-parameters for the d-wave 
and for the AFM phases can be found in Fig. [3] The arrow 
indicates the schematic flow of the parameters when pressure 
is applied, (b) Schematic phase diagram of organic supercon- 
ductors in parameter space of temperature and pressure. 



VMC calculation, the effect is not as big as in the ear- 
lier Gutzwiller-approximation based calculation. In our 
calculation if we consider a fixed and non-zero t'^ instead 
of a variational one, we would in fact get a slope close to 
the one reported by Gan et al. The second difference is 
that our VMC suggests the onset for superconductivity 
to be at [/ = 2.5 for all cases, while Gan et al find this 
phase boundary changing considerable when tuning t' /t. 
We believe that these differences can be attributed to 
the different method and wave- function used. For com- 



parison with the experiments, we plot a schematic phase 
diagram for the layered organic conductors at the right 
panel of the figure. 

In summary, we have presented the results of VMC 
calculations for a recently suggested model for Gossamer 
superconductivity. Our trial wave function has the ingre- 
dient to describe metallic, AFM and SC states. This was 
archived by means of using Jastrow-factors for partial 
Gutzwiller-projection and AFM-weighting. We showed 
that the VMC result is consistent with experiments, and 
supports the previously suggested analytical variational 
calculations qualitatively, as we were able to identify the 
three expected phases, with the help of measurements of 
the order-parameters for AFM and SC. The exact tran- 
sition line between SC- and metallic-phase, and between 
SC- and AFM-phase differs from the one found previ- 
ously. 
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